###########################################################################################
#Replication code for Oliver, Steven, Ryan Jablonski, Justin Hastings. 
#"The Tortuga Disease: The Perverse Effects of Illicit Wealth" Forthcoming International Studies Quarterly (Accepted May 2016)
#Creates Figures 5
#Version 10 June, 2016
#Estimated on R version 3.3.0 (2016-05-03)
###########################################################################################

rm(list = ls())
library(foreign)
library(ggplot2)
library(foreign)
library(MASS)
library(readstata13)
set.seed(123)

setwd("<directory location>")

piracy<-read.dta13("RansomData.dta")
coef_imports<-read.csv("coef_imports.csv")
coef_imports=coef_imports[1,]
cov_imports<-read.csv("cov_imports.csv")
cov_imports=cov_imports[complete.cases(cov_imports),]
coef_imports <- as.vector(as.matrix(coef_imports))
cov_imports=as.matrix(cov_imports)


piracy$yvar<-piracy$log_imports
piracy$xvar<-piracy$log_ransom_amount 

piracy$treat<-piracy$xvar*piracy$region
piracy<-subset(piracy, !is.na(yvar))
piracy<-subset(piracy, !is.na(treat))

#setup a prediction dataset with variables set at their mean
new.data1 = with(piracy, expand.grid(region=mean(region, na.rm=TRUE),
                                     xvar=mean(xvar, na.rm=TRUE),
                                     treat= seq(min(treat), max(treat), by=.1),
                                     log_suezvessels=mean(log_suezvessels, na.rm=TRUE),
                                     log_monthcount=mean(log_monthcount, na.rm=TRUE),
                                     monsoondummy=mean(monsoondummy, na.rm=TRUE),
                                     gdp_growth=mean(gdp_growth, na.rm=TRUE),
                                     ramadan=mean(ramadan, na.rm=TRUE),
                                     livestockban=mean(livestockban, na.rm=TRUE),
                                     pmf=mean(pmf, na.rm=TRUE),
                                     mmmonth_2=mean(mmmonth_2, na.rm=TRUE),
                                     mmmonth_3=mean(mmmonth_3, na.rm=TRUE),
                                     mmmonth_4=mean(mmmonth_4, na.rm=TRUE),
                                     mmmonth_5=mean(mmmonth_5, na.rm=TRUE),
                                     mmmonth_6=mean(mmmonth_6, na.rm=TRUE),
                                     mmmonth_7=mean(mmmonth_7, na.rm=TRUE),
                                     mmmonth_8=mean(mmmonth_8, na.rm=TRUE),
                                     mmmonth_9=mean(mmmonth_9, na.rm=TRUE),
                                     mmmonth_10=mean(mmmonth_10, na.rm=TRUE),
                                     mmmonth_11=mean(mmmonth_10, na.rm=TRUE),
                                     mmmonth_12=mean(mmmonth_12, na.rm=TRUE)
))

#create a second prediction dataset with treat set at zero
new.data0 = new.data1
new.data0$treat=rep(0, length(new.data0$xvar))


nsims=1000
yhat=matrix(nrow=nrow(new.data1), ncol=nsims)
sim.coef=matrix(nrow(piracy), ncol=nsims)

lm.out = lm(yvar ~ region + xvar + treat +  log_suezvessels + log_monthcount + monsoondummy + gdp_growth + ramadan + livestockban + pmf + mmmonth_2 + mmmonth_3 + mmmonth_4 + mmmonth_5 + mmmonth_6 + mmmonth_7 + mmmonth_8 + mmmonth_9 + mmmonth_10 + mmmonth_11 +  mmmonth_12, data=piracy)
thismodel=lm.out
tt <- terms(thismodel)
Terms <- delete.response(tt)

#draw from a multivariate normal distribution for simulating effect
sim.coef=mvrnorm(n=nsims, mu=coef_imports, Sigma=cov_imports)

m.mat1 <- model.matrix(Terms,data=new.data1)
m.mat0 <- model.matrix(Terms,data=new.data0)

#move first column to last (intercept)
m.mat1=cbind(m.mat1, m.mat1[,1])
m.mat1=m.mat1[,-1]
m.mat0=cbind(m.mat0, m.mat0[,1])
m.mat0=m.mat0[,-1]


for(i in 1:nsims){
  m.coef <- sim.coef[i,]
  fit1 <- as.vector(m.mat1 %*% m.coef)
  fit0 <- as.vector(m.mat0 %*% m.coef)
  yhat[,i] = (exp(fit1)-1)-(exp(fit0)-1)
}


#create dataset with simulated predictions
predict.data=data.frame(cbind(upper=apply(yhat, 1, quantile, probs=.95), lower=apply(yhat, 1, quantile, probs=.05), yhat=apply(yhat, 1, mean), treat=seq(min(piracy$treat), max(piracy$treat), by=.1), sd=apply(yhat, 1, sd)))


#transform treatment variable for display
piracy$treat_exp1 = (exp(piracy$treat)-1)/1000000
piracy$treat_exp = piracy$treat_exp1 *(16.9/max(piracy$treat_exp1))

predictplot=ggplot(predict.data, aes(x=treat, y=yhat)) + 
  geom_linerange(ymax=predict.data$upper, ymin=predict.data$lower, colour="grey") +
  geom_line(colour="black",  size=2) + 
  geom_hline(yintercept = 0, size=.8, linetype=2) + 
  geom_rug(data=piracy, aes(x=treat_exp, y=yvar), sides="b", size=.6, position=position_jitter(width=.1, height=0)) +
  theme_bw(base_family="serif", base_size=20) +
  xlab("Ransom Amount") +
  scale_x_continuous(breaks=c(0.01, 16.9), labels=c("MIN\n$0 Million", "MAX\n$23 Million")) +
  theme(
    panel.grid.major=element_blank(),
    panel.grid.minor=element_blank(), 
    panel.border = element_blank(), 
    panel.background = element_blank(),
    axis.line = element_line(colour = "black", size=1),
    plot.margin =  unit(c(1,1,1,1), "cm"),
    axis.title.x = element_text(vjust = 2),
    axis.title.y = element_text(vjust = .25),
    axis.text.x = element_text(size=19,vjust = 0),
    axis.ticks.x = element_line(colour = "black", size=1),
    axis.ticks.length = unit(c(.2), "cm")
  ) 

predictplot = predictplot + ylab("Change in Imports\n(USD Million)") +
  scale_y_continuous(limits = c(-2000000, 8000000), breaks=c(-2000000, 0, 2000000, 4000000, 6000000, 8000000 ), labels=c("-2","0","+2", "+4", "+6", "+8"))
loc="imports_predict.tiff"

predictplot
ggsave(loc, plot = predictplot, width = 7, height = 7, dpi=500)
